Fluctuations of 1/f noise and the low frequency cutoff paradox 
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Recent experiments on blinking quantum dots and weak turbulence in liquid crystals reveal the 
fundamental connection between 1/f noise and power law intermittency. The non-stationarity of the 
process implies that the power spectrum is random - a manifestation of weak ergodicity breaking. 
Here we obtain the universal distribution of the power spectrum, which can be used to identify 
intermittency as the source of the noise. We solve an outstanding paradox on the non integrability 
of 1// noise and the violation of Parseval's theorem. We explain why there is no physical low 
frequency cutoff and therefore cannot be found in experiments. 
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The power spectrum S(f) of a wide variety of phys- 
ical systems exhibits enigmatic 1/ / noise P, 0] at low 
frequencies, 



S(f) 



constant 



7 



where < 7 < 2. 



(1) 



Starting with Bernamont Q , in the context of resistance 
fluctuations of thin films, many models of these widely 
observed phenomena were put forward. Indeed, 1// noise 
is practically universal, ranging from voltages and cur- 
rents in vacuum tubes, diodes and transistors, to annual 
amounts of rainfall, to name only a few examples. A 
closer look at the phenomenon reveals several themes 
which demand further explanation. The first is that 1/f 
noise is not integrable: f_ S(f) df — 00, due to the 
low frequency behavior, when 7 > 1. This violates the 
Parseval theorem from which one may deduce that the 
spectrum of a random process is integrable (see details 
below). So how can we find in laboratory 1/f noise if a 
mathematical theorem forbids it? One simple explana- 
tion would be that the phenomenon has a cutoff at some 
low frequency, namely that below f < fa Eq. (p} is not 
valid. Experimentalists have therefore carefully searched 
for this cutoff, increasing the measurement time as far as 
reasonable: three weeks for noise in MOSFET and 
300 years for weather data [H[ . No cutoff frequency is ob- 
served even after these long measurement times. This is 
one of the outstanding features of 1 // noise. A second old 
controversy, related to the first, is the suggestion of Man- 
delbrot Q that models of 1// noise for 7 > 1 should be 
related to non-stationarity processes, though the nature 
of this non-stationarity is still an open question. Fur- 
ther, experiments find that at least in some cases the 
amplitude of the power spectrum varies among identical 
systems measured at different times, but the shape and 
in particular the value of the exponent 7 is quite consis- 
tent [l], Q • This means that a 1// spectrum is a non self 
averaging observable, at least in some systems. 

While these observations where made long ago, the 
verdict on them is not yet out. However, recent mea- 



surements of blinking quantum dots 0, H| and liquid 
crystals in the electrohydrodynamic convection regime 
Q shed new light on the nature of 1/f noise. Both sys- 
tems, while very different in their nature, reveal a power 
law intermittency route to 1// noise. This means that 
power law waiting times in a micro-state of the system 
are responsible for the observed spectrum. This approach 
was suggested as a fundamental mechanism for 1/f noise 
in the context of intermittency of chaos and turbulence, 
with the work of Manneville [lOj. Subsequently, it has 
been found in many intermittent chaotic systems [lll - 
and has been used success fully as a model for trans- 
port in geological formations [15(. For a quantum dot 
driven by a continuous wave laser, this mechanism means 
that the dot switches from a dark state to a bright state 
where photons are emitted, and that sojourn times in 
both states exhibit power law statistics which is scale 
free [ljl [ljj ■ Waiting times probability density func- 
tions (PDF) in these states follow ip(r) ~ r~( 1+Q ) and 
< a < 1 (for bright state times this behavior is found 
for low laser intensity and low temperature). The dy- 
namics is scale free, because the average sojourn times 
diverge, and we expect weak ergodicity breaking [3, . 
This means that the power spectrum remains a random 
variable even in the long time limit [2l| . 

Here we investigate the non self averaging power spec- 
trum, and show that indeed this observable exhibits large 
but universal fluctuations, while the estimation of 7 is 
rather robust. Our work gives experimentalists a way to 
verify whether a data set exhibiting 1/f noise belongs to 
the intermittency class, and this we believe will help un- 
ravel the origin of an old mystery of statistical physics. 
We also remove the paradox based on Parseval's iden- 
tity, showing that as t — > 00 the integrability remains, 
and that there is no cutoff frequency /q. So experimen- 
tally searching for this "lost" low frequency might be in 
vain. 

Parseval's identity and 1/f noise. We consider a 
measurement of a random signal I(t) in the time in- 
terval (0,t), so that its Fourier transform is It(w) = 



2 



J I(t') exp(— iut')dt'. The power spectrum S t (u) = 
[it (w) i t * (w)]/t, is considered in the long measurement 
time limit. The ensemble average power spectrum is 
(St (w)}. Note that in an experiment with one realization 
of the time series, for example a measurement of the in- 
tensity of a single molecule or a quantum dot, the ensem- 



ble average is not performed, though in experiments one 
introduces smoothing methods which reduce the noise 
level of the reported power spectrum 22] . More impor- 
tantly note that the integral over the power spectrum 
is 



St (w) dui 



dui / dti exp {—iujti 
lo 



•1(h) f 
Jo 



dt 2 exp (iwt 2 ) 1(h) 



2tt 
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/ i 2 (h)dh, 
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where we used a well known identity of the delta func- 
tion f_ dcoexp[— iu)(t\— h)] = 2it5(t2— ti) and St(w) — 
St{—io) by definition. For any bounded process be 
it ergodic or non-ergodic, stationary or non-stationary, 
I 2 {t) < (/max) 2 and hence S t (w) do; < 2^(/ max ) 2 . 
So the integral is finite for a wide class of processes. As a 
consequence, non-integrable 1// noise is strictly prohib- 
ited. The classical way out was to assume a violation of 
the l//-bchaviour in the limit of / — » 0. 

The ensemble-average of Eq. @, under the additional 
assumption that the process reaches a stationary state, 
reads lini{_>. 00 S t (oj)duj = 2tt(I 2 ). If the system is 

ergodic, i.e., I 2 = J* I 2 (t')dt'/t — > (I 2 ), we have for a 
single trajectory lit') 



lim 



St(to)duj = 2ir(I 2 ). 



(3) 



Thus fluctuations of the total area under the power spec- 
trum are an indication for ergodicity breaking. 

Models. For simplicity we consider a two state model. 
The generalization to TV states will be discussed else- 
where. The two states of our model are up or + where 
I(t) = Io and down or — with I(t) = —Jo- The sojourn 
times in these states are independently identically dis- 
tributed random variables with PDFs ip(r). Thus after 
waiting a random time in state up, the particle jumps 
to state down and vice versa. The waiting times PDFs 
have long tails ip(r) « t~( 1+q ) with < a < 1, hence 
the averages of up and down times are infinite. The 
Laplace t —> A transform of these PDFs is for small A: 
■0(A) ~ 1 — (rA) a where f is a scaling constant. This is 
a simple stochastic model of a blinking quantum dot, for 
which typically a = 1/2 though 1/2 < a < 1 was also 
reported. 

More general models where the random process I(t) 
has N internal states with possible different waiting time 
distributions can describe annealed trap models used for 
glass phenomenology or for continuous time random 
walks describing motion of single molecules in live cells 
[23j . Our results for N > 2 are structurally similar to 
those for N — 2 but deserve their own discussion, which 



will be presented in a longer paper. For these models, we 
have derived detailed expressions for the power spectrum 
of the process I(t) and its statistical properties. 

Statement of the main results. For a < 1 the expecta- 
tion value of the spectrum is not constant, but decreases 
with measurement time (S t (oj)) ~ t a ~ 1 a a (ui). Expand- 
ing the ^-independent function a a (to) for small frequen- 
cies id, one finds a typical non-integrable l//-noise 



(S t {w)) * C- 



,2-a ' 



(4) 



In general, the value S't(w) of the spectrum is a fluctu- 
ating quantity even in the t —> 00 limit. The statistical 
behavior of the general class of processes for large t (for 
pairwise disjoint u>i ^ 0) is fully described by 



S t (ui) 



S t (uj v 



(S t (UJ 1 )) , ^^^ , (StK)) 



Y a ■ (£1 , . . . , 



(5) 



where Y a is a random variable of normalized Mittag- 
Leffler distribution with exponent a whose moments are 
(Y"™) = n!r(l + a) n /T(l + no) H. The & are inde- 
pendent exponential random variables with unit mean. 
For a = 1 the Mittag-Leffler random variable becomes 
Y\ = 1, so that the powers St(u>i) of different frequen- 
cies become independent exponentially distributed ran- 
dom variables - a result known for several ergodic random 
processes (25|. In the case of weak ergodicity breaking 
(a < 1), the whole spectrum obtains a common random 
prefactor Y a which shifts the complete observed spec- 
trum. 

Many procedures for the estimation of the spectrum 
from one finite time realization are designed to suppress 
the statistical fluctuations due to the uncorrelated ran- 

13, 



251. These cannot account for the 



dom variables £j 
fluctuations of Y a , common to all estimators of a given 
realization. For these procedures the prefactor affects 
all estimated values for the spectrum. However, being 
a common prefactor, it does not affect the shape of the 
estimated spectrum so that features as l//-noise can be 
detected independently of the realization. 
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Motivation of the results. Whereas an exact derivation 
for the general iV-state model is based on (2(| and will be 
presented in a longer paper, we report here on the 2-state 
model introduced above. To simplify our arguments fur- 
ther, we assume that after a waiting time the next state 
will be chosen randomly (i.e., a + state with probability 
1/2). 

Let Tj be the ith waiting time and Xi = i-^o the 
value taken during this waiting time. We denote by 
Tj = ^21—1 Ti the times at which one waiting time ends. 
If n(t) is the number of completed waiting times up to 
time t, we can approximate by ignoring the waiting time 
in progress at t 



I At exp(zwr)J(r) 
Jo 



n(t) 

i=i 



dj(uj) 



... . . . 1 — exp(iu)Tj) , . 

with dj(ui) = %Xj exp[iuiTj) . (6) 



With this approximation and (XiXj) — ^j^o one obtains 



S t (u) 



n(t) 

- d k (u})di(-uj). 



(7) 



Assuming that n(t) is for large t independent of the wait- 
ing time of a single step, Ti, we get for the ensemble av- 
erage 



(St(u)) ~ ^-(d 1 (Lj)d 1 (-Lj)) 



2 {n(t)) 2 - ip(iw) - ip(-iu) 
~t ~ 2 ' 



(8) 



It has been shown that n(t) ~ Y a t a /(T(l + a)r a ) [22 
Therefore 



(St(w)> 



T 2 f a-1 



ip(iLj) — ip(—ioj) 



T(l + a)T a 
2J^cos(a7r/2) t a ~ l 



(9) 



r(l + a) 



as u> — > 0. 



The last line shows the typical 1/f noise 2l[. It is im- 
portant that the observation limit t — > oo is taken before 
the frequency limit u> — > 0. 

We motivate the main result Eq. ([5]) with help of a ran- 
dom phase approximation. The random phase approxi- 
mation assumes that terms of the form exp(iu;Tj) are just 
random phases and any average over them vanishes. Es- 
pecially, {d jl (i/i) ■ • • dj n (v n )) — if v\Tj 1 H Vv n T ]n / 

for some TjS (the v being ±w). Looking at the second 
moment of Eq. and using (jHJ : 



1 



^(t) 



(s?( w )>~^< ^ d fc Md«(-u>)d P Md,(-w)> 

k,l ,p,q— 1 



f 2 



(7l(t) 2 )(dlMdl(-tj)) 2 



(10) 
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FIG. 1: (color online) a) St(f) plotted for different realiza- 
tions (a = 0.5, t = 10 2 ). Inside each realization one has 
fluctuations following exponential distributions. Different re- 
alizations are shifted with respect to each other due the ran- 
dom prefactor Y a . b) Ensemble average of St(f) plotted for 
different lengths t of the time series. One sees the decay of 
the spectrum (St(f)} =± 0.101 t~ 1/2 f~ 3/2 both in time and 
frequency. The crossover frequency is around f c ~ 0.51/t. 
The simulation perfectly match the theory Eqs. (|4I13[) . 



where we ignored terms with k = l= p = qas there 
are only (n(t)) of them. The factor 2 stems from the fact 
that the sum in the first line of Eq. (fTOj) has contributions 
for k — I, p — q and for k — q, I — p. In contrast to this, 
for the term (St{oJi)St(oJ2)) with cji ^ cu2 this symmetry 
factor will not be present. Following the same steps as 
in Eq. dXOj gives 



(11) 



(S t ( Ul )S t (ua)) c (r a 2 )(^( Wl ))(^(^ 2 )}. 



This shows the equality of the second moments of Eq. ([5]) . 
The equality of the higher moments follows similarly by 
using combinatorial methods to determine these symme- 
try factors. We see that the random number of jumps 
n(t) is responsible for the Mittag-LefHer fluctuations 
while the random phases generate the exponential noise. 
See Supplemental Material at [URL will be inserted by 
publisher] for a more detailed version of this approxi- 
mation. The exact proofs for the general model will be 
published in a longer paper. 

Numerical results. We simulated the two state model 
with Io = 1 for different length of time series and dif- 
ferent a. The waiting times were generated by using a 
uniformly distributed random number < X < 1 and 
setting t = c a X~ x l a . The constant c a was chosen such 
that (n(l)) ~ 10 000. The ensemble consists of 10 000 
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FIG. 2: Distributions of the frequency averaged spectra M 
(see Eq. JUJ). The lines are the analytic probability densities 
of the Mittag-Leffler distributions (t = 10 4 ) 



realizations of the time series. 

In Fig. QJ, we have plotted different realizations of 
St (/)• The stochastic fluctuations inside and between the 
realizations are clearly observable. In Fig. Q~|d the ensem- 
ble average of the power spectrum for different lengths 
is plotted. The 1// spectrum and its decay with obser- 
vation time is clearly visible. Note that at very low fre- 
quencies we find St(w) ~ const independent of frequency 
- an effect we soon explain. 

In a second step we want to check the statistical prop- 
erties described by Eq. ([5]). To isolate the Mittag-Lcfflcr 
fluctuations, we have calculated the spectrum for a fixed 
set of N frequencies u>i and determined the values 



M = 



1 N 

-Y 



(12) 



As the exponential distributions are uncorrelated, they 
average out for sufficiently large N and the value taken 
by M should be distributed as the Mittag-Leffler distri- 
bution Y a . We have compared this for different a values. 
The histogram of M values with the Mittag-Leffler den- 
sity is shown in Fig. [2] for a = 0.2, a = 0.5 and a = 0.8. 
A good agreement with the theory is apparent. 

Removing the non integrability paradox of 1/f noise. 
As mentioned in the introduction, the l/u> 2 ~ a noise is 
non integrable, J Q (S(ui))du> = oo, due to the low fre- 
quency behavior. This in turn violates the simple bound 
we have found. To start understanding this behavior no- 
tice that the random phase approximation breaks down 
when u> = 0, as the phase uT n is clearly non random. 
Hence the distribution of the power spectrum in Eq. ([5]) 



is not valid for u = and this case must be treated sep- 
arately. 

2 

For ui = we have for a single realization St(0) = It 
with the time average / = J* I(t')dt' jt. For ergodic pro- 
cesses the time average I is equal to the ensemble average 
(I). However, for non ergodic processes under investiga- 
tion, the time average / remains a random variable even 
in the infinite time limit 2^, 28 1. For the two state pro- 
cess introduced above, we have I = Iq{T + — T~)/t where 
is the total time spent in state up or down. The value 
of / follows an arcsine like distribution [2(J [28l - [30j . This 
simply means that for a given realization, the system 
will spend most of the time either in state up, or in state 
down, and hence / is random, which would not be the 
case for an ergodic process. 

This has a consequence for the non integrability of the 

2 

power spectrum. As St(0) = I t, the spectrum at zero 
frequency tends to infinity, but for any finite measure- 
ment time it is finite. For the two state model we have 
on average (iSt(O)) = Jq(1 — a)t. So indeed theoretically 
there is a low-frequency cut-off of the divergence of the 
1/f spectrum, and we now define a crossover frequency 
oj c for the transition between the zero frequency limit, 
where arcsine statistics takes control (failure of random 
phase approximation), and higher frequencies where the 
Mittag-Leffler statistics takes control. This frequency is 
defined by merging the two behaviours, 



(St(0)). 



We see that 



_ 2 \l/(2-a) 1 

c/(f)) - 



(13) 



(14) 



The values of (J ) and C can be obtained from mea- 
surements, or from theory, for example for the two state 
model (f) = 1^(1 -a) and C = 2I$cos{an/2)/T(l + a). 
More importantly, we see that the crossover frequency 
depends on the measurement time as 1/t. Although at 
first sight surprising, this is the only way how such a 
cross-over can take place in the absence of a characteris- 
tic time scale for dynamics: measurement time itself sets 
the time scale for crossover. Additionally, 1/t appears as 
the frequency resolution of the discrete Fourier transform 
typically used in spectral analysis. Importantly, experi- 
ments report a lowest frequency at / = 1/t. 

We see that increasing measurement time merely 
stretches the domain of frequency where the 1/f noise 
is observed, which is clearly seen in the numerical simu- 
lations (see Fig. [!})). There is no point in increasing mea- 
surement time in order to identify better the crossover, 
since a time independent crossover frequency does not 
exist. Thus the 1/f noise stretches to the lowest frequen- 
cies compatible with measurement time (of the order of 
1/t). This resolves the non integrability paradox. The 
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amplitude of the power spectrum itself is also decreasing 
in time, in such a way that integrability is maintained. 
Namely 



(S t (w)}dw- (t)tu + I Ct a - x lu 2 - a 



f^((7V- C ) 



l/(a-2) 



(15) 



is indeed finite and time independent. 

Thus we conclude that the power spectrum is inte- 
grable as it should. This seems to indicate the generality 
of our results, since a crossover frequency is only found in 
few experiments. From a different angle, assuming that 
the natural frequency is also the limit of measurement 
uj c ~ 1/t, we must demand the decrease of the amplitude 
of power spectrum with time to maintain integrability 
as required for bounded signals. This together with the 
universal fluctuations of 1// noise (Eq. ([5])) are strong 
fingerprints of power law intcrmittcncy. The tools de- 
veloped here can be tested in a vast number of physical 
systems. 
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Supplementary material to: 
Fluctuations of 1/ f noise and the low frequency cutoff paradox 

Markus Niemann, Holger Kantz, and Eli Barkai 



This supplementary material contains a more detailed derivation of the main results. We will still work with the 
^sl \ approximations made in the main text - an exact proof will be published in a longer paper. 
Using the same notations as in the main text we see 

dj(ui) = / dr Xj cxp(icjr) 

(SI) 



O 

(N 



<Z3 



O 



X 



T, 



i ■ m \ 1 — exp(iwr 7 ) 



,t n(t) Tj+1 t 

/ dr exp(zwr)/(r) = 2^ / dr exp(iwr) + / dr Xj exp(iwr) 



53: 

= exp(iwTj-) 

With this, we get for the Fourier transform of the signal by neglecting the last waiting time in progress at t: 
■ (S2) 

g n(t) V ' 

where n(i) denotes the number of finished waiting times up to time t also called the number of renewals in (0,i). 
With this approximation we get for the spectrum 

1 Mt) 

S t {uj) ~ - d k (u)di(-u). (S3) 
fcj=i 



We now assume that the Ti are independent from the tj and in turn the average of n(t) can be performed independently 
from the average over the Tj. Using (xiXj) — ^ij^o one obtains (where a distinction is necessary, we write the variables 
whose average is taken as an index) 



n(t) 

<S t H>~ 7 (^ d*( w )d,(-a;)> 



1 / / \ / / * frp T ul" exp(io;r fc ) 1 - exp(-icjr ; ) 



C<i '■ _ T 2 1 /y' ( 2 ~ exp(?urr fc ) - exp(-zo;rfc)) Ti \ ^ S4 ' ) 

\fc=l / n(t) 



S 1 „ (»(*)> 



2 {n(t)} 2 — tj){iw) — 

n : 2 

(di(w)di(-w)). 



10 t 



t 

It has been shown that n(t) ~ Y^/O^l + a)r Q ) from which we get (n(t)) ~ r/(r(l + a)r a ) 

a f;-' (S5) 

1(1 + a)r ur 



Here Y^, is a random variable of Mittag-Leffier distribution with unit mean which is defined by its moments (Y™) — 
n!F(l + a) 2 /r(l + na). Using a random variable L a of totally asymmetric Levy stable distribution with exponent a 
and distribution i a (£) 

<exp(~AL Q )) = / d£ exp(-A£)U0 = exp(-A Q ), (S6) 
Jo 
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a representation of Y a can be obtained by 

Y a = T(l + a)L- a 
which is equivalent of Y a having the probability density function 0, 0| 



r 1 / Q (i + a) 



(S7) 



(S8) 



Plots of y a (0 are used in Fig. 2 for various a. 

For small A we have the expansion in the Laplace variable tp(X) — 1 — (rX) a . By analytic continuation (putting the 
branch cut to the negative real axis such that it is not crossed) we have 



exp 



i-sgn(uj) 



(*«r = (M. 

= \oj\ a (cos (~a) + i sgn(w) sin (~c 



such that for small w 



and 



1 — ijj(iuj) ~ r a |cj| Q (cos (^ a ) + * s g n ( w ) sm (t^c 



2^cos(a7r/2) i c 



r(l + a) 



as w — >• 0. 



w ■ 



(S9) 



(S10) 



(Sll) 



It is important, that the observation limit t — > oo is taken before the frequency limit u> — > 0. 

The Eq. (j5j) is shown by checking the equality of the moments on both sides. Assuming that we have some natural 
numbers pi, . . . ,p q (P = p\ + ■ ■ ■ +p q ), we need to show 



Stiujj) 



nv(5 t K-)> 



\i=i 



(S12) 



Here, we denote by = an equality which needs to be proven. We have used the fact that the £j are independent 
random variables which are exponentially distributed with unit mean, i.e., (£") = n\. Especially, we have to show for 
the second moments (wi ^ LO2) 



St(u>) 



= 2(i; 2 ) 



and 



5 t (a;i) 5 t (^2) \ ^ , y2 \ 



(S13) 
(S14) 



We motivate this with help of a random phase approximation. The random phase approximation assumes 
that terms of the form exp(iwTj) are just random phases and any average over them vanishes. Especially, the 
phase factor of dj x (y{) ■ ■ ■ dj n (u n ) is exp (i[viTj 1 + • • • + v n TjJ\) (the take values of ±w). Therefore, we assume 
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(dj 1 (vi) ■ ■ ■ dj n (v n )) = if v\Tj 1 + • • • + v n T.j n ^ for some TjS. Looking at 

n(t) 

<5 f 2 H)^-( ]T d k (Lj)dl(-u)d p (u)d q (-u)) 
p,q=l 
n(t) 

t 2 



1 / ™ (t) \ 

72 \ H (4Mdz(-w)d p (w)d g (-w)) Xi , Ti , ri \ 

\kd,p,q=l I , 



(t) (S15) 

1 / / \ / 1 — exp(iwrfe) 1 — exp(— ilotA 1 — exp(iwr„) 1 — exp(— ituT q ) 
72 \ 1^ \XkXlX P Xq)xi { 

t A \ ' \ LJ UJ UJ UJ 



(exp(«o;(r fc - Ti + T p - T q ))) Ti 



i(t) 



In our approximation, the non vanishing terms in the sum must have the property u){T k — Jj + T p — T q ) = for any 
realization. This is possible either for k = I and p = q, or for k = q and p = I. In other words, the random phase 
approximation boils down to 

(exp(«o;(r fc - Ti + T p - T q ))) T( ~ S k iS pq + 8 kq 8 p i. (S16) 
Plugging this back into Eq. (|S15[) and using (fS4|) gives 

2 / n(t) \ 1 /™ (t) \ 

<s t 2 H) ^ s ( E (4(^)4(-^))(d P (^)d P (-w)) ) + ^ ( ^(4(w)4(-w)4M4(-w)) ) 
\fc,p=i / \fe=i / 



2 , ,. w 1 

|(n(t) 2 )(rf 1 (u;)d 1 (-a;)} 2 



~<n(t)(n(t) - l))(d 1 Md 1 (-o;)) 2 + ^<n(t))(di(w)di(-w)di(w)di(-w)) ( S17 ) 



^2(r Q 2 )(5 t M) 2 . 

where we ignored terms of the order (n(t)) as for large t the leading order is (n(i) 2 ). This shows Eq. (|S13I) . 
In contrast to this, for the term (St(uji)St(uj2)) with ui\ ^ cj 2 we get 

n(t) 

(S t (uii)S t {u) 2 )) ~ -j( 2J d/c(wi)di(- w i) rf p( w 2)d 9 (-^2)) 

k,l.p,q— 1 



1 / n(t) \ 

72 \ X! ( rf fc(wi)d;(-Wi)dp(w 2 )dg(-W2)) Xi ,T i ,r i ) 

\fc,J,p,g=l / „ (4) 

1 / , , / 1 - exp(iwiTA;) 1 - exp(-«wir ; ) 1 - cxp(iw 2 T„) 1 - exp(-zw2r„) 

72 ( V (XkXiXpX q )xt ( " 

V \ ' \ wi wi w 2 w 2 

k,l, P ,q— 1 



(exp(iw 1 (T fe - T ( ) + iw 2 (T p - T ? ))) 



(S18) 

This time, the non vanishing terms in the sum must have the property uj\{Tk — TJ) +w 2 (T p — T q ) — which is possible 
only for k — I and p = q: 

(exp(ia;i(T fe - Ti) + iui 2 (T p - T q ))) Ti ~ 5 (S19) 
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We get with this approximation 

1 / «(*) V 1 /»(*) \ 

(S t (u>i)S t (uj2)) - Ta ( ^ (<4(wi)t4(-Wi))(d p (w 2 )dp(-W2)) ) + 72 ( ^( d fc( w i) d fc(- w i) rf fc(w 2 )d fc (-cj2)) / 
\fc,p=i / \fc=i / 

~ ^(n(t)(n(t) - l)><di(wi)d 1 (-wi)))<di(w 2 )di(-W2)> + ^(n(t))(di(wi)di(-wi)di(w 2 )di(-W2)) 

^(y Q 2 )(5 4 ( Wl ))(5 t (^)). 

(S20) 

This shows Eq. (|S14[) and therefore with Eq. (|S 13[) the equality of the second moments of Eq. ([5]) . Similar calculations 
can be performed in general for higher order moments giving rise to Eq. ([3]). 
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